%% EXERCISE 22: SECOND ORDER PHASE FIELD MODEL - students's version

clear all
%close all
%clc

format short
addpath('ex22_toolbox_local')


% Nodal Spacing options
optNS.dim = 2;
optNS.knn = 0;

%% penalty parameter
p_val = 10000;

%% LME options (Machine precision)
optLME.dim   = 2;
optLME.verb  = 0; % 0:off    1:on
optLME.knn   = 0;
optLME.grad    = 1;           % Computation of the Gradient 0:OFF 1:ON
optLME.hess    = 0;           % Computation of the Hessian  0:OFF 1:ON
optLME.TolNR   = 1.e-12;      % Newton-Raphson Tolerance
optLME.Tol0    = 1.e-8;       % Target Zero Tolerance

gamma = [4.8]; %LME gamma aspect parameter [Arroyo, IJNME 2006]
%gamma = [0.8 1.8 4.8]; <-  For exercice 2.2.3

%% Integration options
%Numerical integration, cubature order
optGL.orderGL = 2;  %order 2(3 gPts), 3(4), 4(6), 5(7), 6(12), 10(25)
optGL.quality = 0.01;

%% Geometrical Parameters
Lx    = 1.0;    % length
Ly    = 1.0;
Ndof  = 21;

tic
disp('------------------------------------------------')
fprintf(1,'Grid size: %2dx%2d\n',Ndof,Ndof);

% Node Points
Nx      = Ndof;
[x_nodes x_nodes0 id_bd] = SquareNodes(0, Nx, Lx);
nPts0   = length(x_nodes0);
nPts    = length(x_nodes);
Ndof    = Nx+2;
h   = Lx/(Nx-1);
%h   = Lx/(Nx-1);
h_n = ones(nPts,1)*h;

optLME.spacing= h;
optNS.spacing = h;
optLME.h_n    = h_n;

%% Phase-field model parameter
ele_n = (0.84)*h; %[0.5 1 2 4]*h; -> 0.84 Best
G_n   = zeros(length(gamma),length(ele_n));

% GAUSS_LEGENDRE
[x_samples w_samples] = MakeGLSamples2D(x_nodes0, optGL);

%% LME: The displacement field is computed
for g=1:length(gamma)
  optLME.gamma = gamma(g);
  t=cputime;
  fprintf(1,'\tgamma=%4.2f | order=%2d\n',optLME.gamma,optGL.orderGL);
  
  % Nodes thermalization, nodes and samples adjacency structures are computed
  [beta_n range_n] = NodalThermalization(x_nodes, optLME);
  optLME.beta   = optLME.gamma/(h*h);
  optLME.beta_n = beta_n;
  optLME.range_n= range_n;
  
  % The basis functions are computed
  % adjacency structure with the nearest neighbors nodes to each sample point
  % nodal shape parameter
  s_near = SamplesAdjacency(x_nodes,x_samples,range_n);
  
  % Local-max entropy basis functions computation
  optLME.s_near = s_near;
  outLME = wrapper_lme(x_nodes,x_samples,optLME);
  
  p_samp  = outLME.p_samp;
  dp_samp = outLME.dp_samp;
  fprintf(1,'cputime LME   : %4.3f\n', cputime-t);
    
  for le =1:length(ele_n)
    ele = ele_n(le);
    % SOLVE SYSTEM
    % ASSEMBLY: A = M + ele^2 * K.
    %The mass matrix M and the stiffness matrix K are computed
    t = cputime;
    optSystem.optLME  = optLME;
    optSystem.pdeOrder= 2;
    optSystem.ele     = ele;
    optSystem.p_val   = p_val;
    optSystem.nPts    = nPts;
    optSystem.spacing = h;
    optSystem.s_near  = s_near;
    optSystem.p_samp  = p_samp;
    optSystem.dp_samp = dp_samp;
    optSystem.w_samp  = w_samples;
    
    %%
    v_lme = ex22_SolveSystem(id_bd,x_nodes,optSystem);
    fprintf(1,'cputime System: %4.3f\n', cputime-t);
    
    % Energy
    t = cputime;
    G = 0;
    for k =1:size(x_samples,1)
      k_near = s_near{k};
      u_k    = v_lme(k_near);
      p_k    = p_samp{k};
      dp_k   = dp_samp{k};
      
      %    Exercise 2.2.1
      u_h   = p_k'*u_k;
      du_xh = dp_k(:,1)'*u_k;
      du_yh = dp_k(:,2)'*u_k;
      
      %    Exercise 2.2.1
      G = G + 1/(2*ele) * (u_h.^2 + ele^2*(du_xh.^2 + du_yh.^2)) * w_samples(k);
    end
    %    Exercise 2.2.1
    G_n(g,le) = G;
    fprintf(1,'cputime Gamma : %4.3f\n', cputime-t);
  end
  disp('------------');
end
disp('---------------------------------------------------');

%% Plot crack length
G_0   = 0.5;
xdata = ele_n/h;

G_sol = G_0./G_n;

figure(1)
plot(xdata,G_sol(1,:),'ro-');
xlabel('ele/h','fontsize',24,'fontweight','b')
ylabel('\Gamma/\Gamma_L','fontsize',24,'fontweight','b')
legend(strcat('\gamma = ',num2str(gamma(1))));
title(strcat('\Gamma/\Gamma_L - 2D - 2nd order PDE - h=',num2str(h)),'fontsize',24,'fontweight','b')
ylim([0,1.2])

%% Plot numerical solution
figure(3);clf
nSY = 51;
nSX = 51;
plotLME2D(x_nodes,v_lme(1:nPts), nSX, nSY, Lx, Ly, [0 0], optLME);
view([-17 4])
title(strcat('LME  \gamma = ',num2str(gamma)),'fontsize',18,'fontweight','b')
xlabel('X', 'fontsize',16,'fontweight','b')
ylabel('Y', 'fontsize',16,'fontweight','b')
zlabel('\delta_z','fontsize',16,'fontweight','b')
fprintf(1,'cputime Plot: %4.3f\n', cputime-t);